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Abstract. - We compare the harmonic and anharmonic properties of the vibrational modes in 
3-dimensional jammed packings of frictionless spheres interacting via repulsive, finite range poten- 
tials. A crossover frequency is apparent in the density of states, the diffusivity and the participation 
ratio of the modes. At this frequency, which shifts to zero at the jamming threshold, the vibra- 
tional modes have a very small participation ratio implying that the modes are quasi-localized. 
The most anharmonic modes occur at low frequency which is opposite to what is normally found 
in crystals. The lowest frequency modes have the strongest response to the pressure and the lowest 
energy barriers to mechanical failure. 
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Introduction. — A starting point for understanding 
dny solid is a calculation of its harmonic vibrational exci- 
tations. However, many important properties require an 
linderstanding of anharmonic effects as well. For exam- 
ple, in a crystal, heat transport and thermal expansion 
are governed by the anharmonic coupling of the harmonic 
modes [1]. Anharmonic effects become even more essen- 
tial when one considers how a solid, be it crystalline or 
amorphous, disintegrates and loses rigidity. Such behavior 
requires the system to depart from the harmonic regime 
as it continually moves over potential-energy barriers to 
explore different configurations. In this paper we discuss 
the anharmonic properties and energy barriers associated 
with the normal modes of an amorphous solid compressed 
above the jamming threshold. 

At the jamming transition, the system is precariously 
perched between a solid and a liquid. The vibrations 
calculated for such a solid composed of soft frictionless 
spheres interacting via repulsive forces possess decidedly 
unusual properties [2] . Rather than having a Debye spec- 
trum in which the density of normal modes, D{uj), at low 
frequency, oj, varies as: D{ui) cx u!'^~^ in d dimensions, the 
density of states at the transition packing fraction, (j)c, has 
a plateau that extends all the way down to zero frequency. 



Thus, there are many more low- frequency modes than can 
be accounted for simply by plane- wave excitations. Upon 
compression to a packing fraction, (j) > (j)c, the plateau 
persists but only down to a cutoff frequency, uj*. Similar 
behavior is found for ellipsoid packings [3] and systems 
with friction [4] . It is known that various geometrical fea- 
tures of the modes depend strongly on frequency [5]. 

Here we show that as the frequency is lowered, the 
modes become progressively more heterogeneous with a 
lower-than-average mode coordination number. Near uj*, 
the modes begin to have high displacement amplitudes in 
small, confined regions of space. Upon compression or ap- 
plication of shear stress, such modes can eventually desta- 
bilize the system when their frequency passes through 
zero. These lowest-frequency modes also have the lowest 
energy barriers to particle rearrangements; even at low 
temperature, they can be sufficiently excited to force the 
system into different configurations. On approaching the 
unjamming transition, the barriers between nearby config- 
urations shrink to zero so that anharmonic effects become 
more pronounced. Thus not only is the jamming transition 
marked by critical behavior in the harmonic properties of 
the solid [6] but also by diverging anharmonic effects. 
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Methods. — We analyze the same system studied pre- 
viously [2,6-8]. Our jammed packings consist of N fric- 
tionless spheres of mass to in a three-dimensional box with 
periodic boundary conditions with 500 < N < 10000. To 
avoid crystallization, wc use a binary mixture of parti- 
cles with diameters ctl and as with ratio aL/crs ~ 1-4. 
The particles interact via a repulsive harmonic potential 
V{rij) = e (1 — Tij /(Tij)^ /2 when the distance between the 



centers of particles i and j , r.y , 



is smaller than the sum of 



their radii, aij, and otherwise. The jammed solids were 
generated through conjugate gradient energy minimiza- 
tion [9]. We tuned the volume fraction and successively 
minimized the energy until a jammed solid at A(/) = — 0^, 
was obtained. We use units where e = 1, CT5 = 1, and 
TO = 1. 

For comparison, we also constructed an unstressed 
counterpart for each configuration studied. Such systems 
offer special insights. To obtain unstressed systems, we 
replaced the interaction between each pair of particles by 
an unstrctchcd spring which has the same spring constant, 
V" ^ as in the original system. Thus, the pressure is au- 
tomatically zero. In comparison to the original "stressed" 
system, the modes of the unstressed one are moved to 
higher frequencies [8,10,11]. 

We diagonalize the Hessian matrix of the jammed solids 
using ARPACK [12] to obtain the frequency u)„ and po- 
larization vector e„^i of each particle i in mode n. From 
these quantities, we calculate the density of states, D{uj), 
and the energy diffusivity of each mode, d{ujn) [7,8], which 
quantifies how rapidly a wavepacket narrowly peaked at a 
frequency cOn spreads out in time. We also calculate geo- 
metrical properties of each mode. The participation ratio 
p{uJn) quantifies how extensive a mode is and ranges from 
(locahzed) to 1 (extended): 



(1) 



We also compute the mode-average coordination number, 
z([j„), which is the effective coordination number weighted 
by the polarization vectors in each mode: 



(2) 



where is the number of interacting neighbors for particle 
i. 

We quantify the anharmonicity of each mode in vari- 
ous ways. We compute the Griineisen parameter, 7(w„), 
which measures the response of the mode frequency to 
compression: 

d(lnw„) 



d{lnq 



(3) 



In jammed solids, both the pressure P and coordination 
number z vary with compression and affect the frequency 
of the modes [6, 10, 11]. Therefore, Eq. (3) is rewritten as 



dP 



dP 
90 



UJ„ 



dz 



(4) 



Near there is always a small fraction ^ of particles 
that are rattlers. The rattlers may become network parti- 
cles upon compression and contribute to 7^: ^zif^n) = 
:^[(l-0(^)i+^(^)2]plf' -here {^)^ and 
(%^)2 ^'"^ change of ujn with respect to z without and 
with the change of the number of rattlers, respectively. In 
our simulations, we measured 7p in the linear response 
regime without changing the number of particle interac- 
tions. We measured 7^ by adding unstretched springs with 
the spring constant V" between network neighbor parti- 
cles or the rattler and its neighbors that do not interact 
but tend to touch under compression, without disturbing 
the pressure. 

We also compute u„iax{^n) which is how far the nth nor- 
mal mode can be displaced (leaving all the other modes 
uncxcited) before the system crosses an energy barrier of 
height Vmax{^7i) and lands in a different configuration. 
To measure Umaxi^n) and Vmax{^n)i we first displaced 
the jammed state at Rq = (ri, r2, r7v)o in the configu- 
rational space along e„ to R(u,a;„) = Rq -|-ite„, and then 
quenched the perturbed system to the local energy min- 
imum. When u < Umaxi'^n), the local energy minimum 
that is obtained is identical to Ro, but when u > Umaxi'^n) 
the local energy minimum changes. Thus, Umax{^7i) sets 
the maximum displacement below which R remains in 
the same energy basin in configurational space. The en- 
ergy barrier of the basin of attraction along mode n is: 

Vmax{(^n) = V {Ii{Umax,U}n)) ~ V {Hq) , whcrC F(R) is tllC 

potential energy per particle at R. Because the jammed 
solid configurations are amorphous, there is no symmetry 
and Umax and Vmax will be different for the -f and — direc- 
tions of each mode. We define Umaxi'^n) and Vmaxi'^n) as 
the minimum values over these two directions. It is impor- 
tant to realize that these quantities do not fully describe 
the extent of the energy basin since linear combinations of 
the different normal modes will give different magnitudes 
for the maximum displacements and energy barriers. In 
general, the basins have a very complex boundary [13]. 

Nature of the vibrational modes: harmonic 
properties. — In Fig. 1, we compare several proper- 
ties of the modes as a function of frequency for both an 
unstressed (left column) and a stressed (right column) con- 
figuration at the relatively large compression of A(/) = 0.1. 
We first examine the harmonic properties: the density of 
states D{ijj) [2], the energy diffusivity d{ui) [7,8], the par- 
ticipation ratio p{u)) [5] and the mode-average coordina- 
tion number z{ijj). 

In the unstressed regime, there are three regimes as 
observed previously for stressed systems [5]: (i) At low 
frequencies there is a regime of low scattering where the 
modes are weakly-scattered plane waves. These modes 
appear as peaks in the density of states, the diffusivity 
and the participation ratio, since they are quantized by 
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Fig. 1: Spectra of (a) density of states D{ui), (b) energy dif- 
fusivity d{uj), (c) participation ratio p{to), (d) mode average 
coordination number z{u}), (e) pressure contribution of the 
Griineisen parameter 7p(tj), (f) coordination number contribu- 
tion of the Griineisen parameter 72(1^), and (g) total Griineisen 
parameter 7(0)) = 7p(a;) +^z{u}) for an unstressed system (left 
column) and a stressed system (right column) at A'^ = 2000 
and A<j) = 0.1. The dashed line in the left column is the guide 
to eye. The solid horizontal line in (gl) shows 7 = 0. Note the 
large difference in scale of the vertical axes in the left and right 
columns for the three types of Griineisen parameter. 



the finite size of our system, (ii) At intermediate frequen- 
cies, above a characteristic frequency (denoted by a ver- 
tical dashed line in the left column of Fig. 1) that scales 
as u!* oc A0°-^, there is strong scattering of the plane 
waves [7,8] and the modes are anomalous and built out of 
the floppy modes at the jamming transition at (j) = (j)c [10]. 
In this regime, the density of states is high and the energy 
diffusivity is small and nearly constant, (iii) At high fre- 
quencies, the modes are localized. 

It is clear from the left column of Fig. 1 that for un- 
stressed systems, signatures of the crossover at uj* appear 



in D{uj), d{uj) and p(w). Previously, it was shown that the 
rise in the density of states at w* corresponds to a boson 
peak [2] and that the leveling off of the energy diffusivity 
at Lo* corresponds to an loffe-Regcl crossover in the energy 
transport [7,8]. We also plot the participation ratio for 
every mode. For the unstressed systems, p{uj„) has a pro- 
nomiced dip at lo* indicating that the modes in this region 
are more localized than elsewhere in the spectrum (ex- 
cept for the Anderson- localized modes at high frequency.) 
Above LO* , the modes are extended with a roughly con- 
stant participation ratio p{lo) « 0.5. The drop in the 
participation ratio near lo* is perhaps not surprising if one 
considers that lo* marks a dramatic drop in the density of 
vibrational modes, which can lead to localization. In this 
case, plane waves below lo* hybridize with these modes, 
so that the resulting vibrations are not truly localized but 
rather are highly confined in small spatial regions with a 
small plane-wave background elsewhere. We will refer to 
these modes as "quasi-localized" modes. In terms of a 
scattering picture, the existence of quasi-localized modes 
at a particular frequency in the soft potential model [14, 15] 
has been shown to lead to strong scattering of vibrations 
at higher frequencies. This is consistent with our find- 
ing quasi-localized modes at the loffe-Regcl crossover fre- 
quency. Fig. Id shows that these quasi-localized modes 
are concentrated in regions where the number of interact- 
ing neighbors is lower than average. The mode-average 
coordination number, z{lo), has a dip at lo* , implying that 
large particle displacements in the quasi-localized modes 
predominantly exist in low-coordination regions. 

The characteristic behaviors described above for D{lo), 
d{Lo), p{uj) and z{lo) at lo* are observed at all com- 
pressions in which lo* falls within the observable range 
for our simulations. In our finite-sized systems, 
plane waves exist at the quantized frequencies LOp^ = 

'^'^(^T,L^Jrt^i'+^^^^'+^^^/ L, which shift down with increas- 
ing system size L cx N^^^ , where ct,l is the speed of trans- 
verse or longitudinal sound, n^, Uy, Uz = 0, ±1, ±2, and 
L is the length of the system. The system size therefore 
limits the observable range of frequencies. At smaller com- 
pressions, where lo* lies below our accessible range, we see 
only regimes (ii) and (iii), corresponding to the plateau in 
the diffusivity and the localized modes at high frequency. 

The modes in a stressed system, as compared to its un- 
stressed counterpart, are pushed to lower frequency since 
repulsive forces tend to lead to buckling instabilities of the 
force network. Comparing the right- and left-hand pan- 
els of Fig. Ic, we find that the effect of pressure on the 
low-participation ratio (quasi-localized) modes is to push 
them towards lo = 0. As a result, regime (ii) is pushed 
down in frequency, so that regime (i) is presumably no 
longer visible. 

The participation ratio data are consistent with earlier 
studies of the frequency-averaged participation ratio [5]. 
We show here the value of p{ujn) for every mode to dis- 
play fine structure that is not easily observed when p(u;„) 
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Fig. 2: (Colour on-line) Participation ratio, p{u)), of stressed 
systems at = 10000 and A(/> = IQ-" (black circles), 10"^ (red 
squares), and 10^^ (blue triangles). Note that only the lowest 
4000 modes are shown in each case. 

is averaged within frequency bins. In Fig. 2, we show the 
participation ratio for the stressed systems at three com- 
pressions. In all cases, p{ujn) decreases dramatically and 
appears to reach the same low value independent of the 
packing fraction. The only difference is that for higher 
values of A0, the anomalous modes do not persist to as 
low a frequency as they do for systems closer to the jam- 
ming threshold. Thus, the excess low-frequency modes, 
while not plane- waves, are also not all uniformly extended. 
These results suggest that as the packing fraction of the 
system is decreased towards the unjamming transition, 
quasi-localized modes always exist but shift to lower and 
lower frequency. 

Because we are dealing with finite-size samples, it is 
difficult to ascertain whether the participation ratio can 
become arbitrarily small and approach zero in the thermo- 
dynamic limit. It is also difficult to determine the prop- 
erties of the low-frequency modes in large systems where 
plane-wave modes are plentiful at low frequency. Since 
there is no symmetry to forbid it, the plane waves will 
hybridize with the excess modes to produce some small 
plane-wave character at large distances. In Fig. 3, we 
show the participation ratio at different system sizes N . 
At w w LOpim there are plane- wave peaks superimposed on 
a smooth background. The peaks shift downwards and 
increase with increasing N , but the background does not. 
This A^-independent background clearly shows the trend of 
decreasing participation ratio with decreasing frequency, 
suggesting that in the large N limit, the participation ra- 
tio at Lo* should approach arbitrarily small values at low 
frequencies. 

Nature of the vibrational modes: anharmonic 
properties. — Fig. 1 also shows one measure of the 
anharmonic behavior of the modes, the Griineisen param- 
eter, versus frequency. We separately plot the pressure 




Fig. 3: (Colour on-line) Participation ratio p{uj) of stressed 
systems at = 0.5 and A'' = 10000 (black solid), 2000 (red 
dashed), 1000 (blue dot-dashed), and 500 (green dotted). Data 
were averaged over configurations and frequency bins. 



and coordination-number contributions to the Griineisen 
parameter 7p(w) and 7z(w) along with the total ^{lo) = 
7p((x') -|- ^z{'^) for both unstressed and stressed systems. 

It has been shown previously that the frequency of dif- 
fusive modes in regime (ii) is determined by the extra co- 
ordination number above isostaticity Az and pressure P: 
up' ~ A(Az)^ — BP, where A and B are positive con- 
stants [10, 11]. We can thus derive the frequency depen- 
dence of the Griineisen parameters: 7p(w) ^ — tj"^, and 
7z(a;) ~ . These scalings are verified in Fig. 4. The 
different frequency dependences of 7p and lead to a 
crossover frequency lo^ at which 7 = 7p -I- 7z = 0. For 
Lo > u}^ the modes move upwards in frequency upon com- 
pression, while for < tj^ the modes shift downwards. 

The inset to Fig. 4 shows estimated at various vol- 

1/2 

ume fractions for stressed systems: ~ (^0) ■ This 
is the same scaling as was found for w*, the onset of ex- 
cess modes (the boson peak) in the density of states [2], 
and for ti^/p, the onset of the low-diffusivity plateau [7,8]. 
Fig. 1 shows that below w-,,, the Griineisen parameter be- 
comes enormous and negative in stressed systems. Thus 
the low-frequency modes are highly anharmonic and un- 
stable. This is completely different from the situation for 
crystals, where the Griineisen parameter is small and inde- 
pendent of frequency: 7(w) = ^^^pn^^^^ ^ dimensions, 
where c(0) is the speed of sound. For a 3-dimensional crys- 
tal of spheres interacting with finite-ranged harmonic re- 
pulsions, 7(0;) = —1/6. Finally we note that the localized 
modes at high frequencies do not show large anharmonic- 
ity, in agreement with the observations in amorphous sili- 
con [17]. 

Fig. 5 shows another measure of the anharmonicity: the 
maximum displacement of Umax and energy barrier Vraax 
along each mode above which the perturbed system falls 
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Fig. 4: (Colour On-liue) Griineisen parameters from the pres- 
sure contribution 7p (circles) and the coordination number 
contribution 7^ (triangles) for stressed systems at A'' = 2000 
and A(j) = 10~*. The red solid and blue dashed lines have a 
slope of —1 and —2, respectively. Data were averaged over con- 
figurations and frequency bins. The inset shows u-y (black solid 
circles), the crossover frequencies at which 7 = 7p-|-7z = 0, as 
a function of Acp for the stressed systems. The red solid line 
has a slope of 1/2. 




out of the original energy basin in configurational space. 
In Fig. 5, results are shown for seven different basins. They 
all follow the same trend in frequency but with signifi- 
cantly different magnitudes. We have also measured these 
quantities at different packing fractions and found similar 
trends in the frequency dependence. Note the significant 
variation in the overall magnitude of Umax from configu- 
ration to configuration. 

Fig. 5a shows that the maximum amplitude of a mode, 
Umax, or the extent of the basin in the direction of the 
mode, increases with increasing lu. This trend is consistent 
with the existence of a large negative Griineisen parameter 
at low frequencies below cj-y. However, it is the opposite 
to what we find in crystals, where the lowest frequency 
plane waves have much larger values of Umax than those at 
higher frequencies. The energy barrier, Vmax, is shown in 
Fig. 5b. Note that the barrier height shows a very strong 
frequency dependence, with the lowest-frequency modes 
having the lowest energy barriers. Thus, at low temper- 
atures, the low-frequency quasi-localized modes are the 
ones that will most easily be driven unstable by thermal 
energy. 

Conclusions. — We studied the harmonic and an- 
harmonic properties of the vibrational modes in jammed 
packings of frictionless spheres. We find that the partici- 
pation ratio of the modes drops with decreasing frequency 
so that modes become highly resonant, or quasi-localized, 
at a characteristic frequency. This frequency decreases as 
the compression is decreased towards the jamming thresh- 
old, so that just above the threshold, the modes appear 
to exist but at nearly zero frequency. Their participation 



Fig. 5: (Colour on-line) Spectra of (a) the maximum displace- 
ment Umax and (b) energy barrier Umax along each mode, n, 
before falling into another basin of attraction for stressed sys- 
tems at A'^ = 1000 and Afji = O.f. Seven configurations are 
shown to demonstrate that the trend in frequency is simi- 
lar in all cases but that there is significant variation in the 
magnitudes of Umax and Umax between different realizations of 
jammed systems. Data were averaged over frequency bins. 



ratio remains low when the system size is increased. In 
some ways, these modes resemble the localized modes at 
the high-frequency end of the spectrum, in that both types 
of modes have small participation ratio and are localized 
in regions with fewer constraints than average. However, 
the quasi-localized modes are significantly different from 
the localized ones in that they are the most anharmonic 
modes in the system, with the largest negative Griineisen 
parameters and the lowest energy barriers for crossing into 
different basins in the energy landscape. 

Similar quasi-localized modes at low frequencies have 
been observed in simulations of amorphous silicon [17-f 9], 
Lennard- Jones glasses [16,20], and soft repulsive disks [21] 
and in experiments on amorphous polymers [22] . In silica, 
the modes also occur near the boson peak frequency [f9], 
the high-displacement regions in these modes are prefer- 
entially situated on under-coordinated atoms [18], and 
the modes themselves have strong anharmonic correc- 
tions [17, 19], consistent with our findings. Our results 
suggest that the existence of highly anharmonic, low- 
frequency, quasi-localized modes near the boson peak fre- 
quency is not just a feature of these specific systems. 
Rather, it is a generic feature not only of disordered net- 
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work glasses (of which our unstressed system is a simple 
example) but also of glasses stabilized by repulsive interac- 
tions, as exemplified by the stressed sphere packings stud- 
ied here. Moreover, since these modes appear to approach 
zero frequency at the isostatic unjamming transition of 
these disordered systems, they can be viewed as originat- 
ing from the soft modes there. 

The strong anharmonicity and low energy barriers asso- 
ciated with the low-frequency quasi-localized modes imply 
that when the system is compressed further or sheared, it 
is these modes that eventually go unstable and give rise 
to initially localized rearrangements in the sample. Thus, 
the deformations in the sample should be confined to small 
regions, as postulated in the picture of shear transfor- 
mation zones [23]. Recent numerical studies on model 
glasses suggest that low-frequency quasi-localized modes 
are correlated with irreversible particle rearrangements at 
nonzero temperature [20] . The low energy barriers that wc 
have found to be associated with the quasi-localized modes 
imply that they should have important consequences for 
system stability even at low temperatures. This suggests 
that further work should be done to examine the rela- 
tion of these modes to the physics of two-level systems, 
which are believed to be the dominant excitations in low- 
temperature glasses. 

Finally, as noted above, the energy barriers of soft- 
sphere systems decrease to zero as the unjamming tran- 
sition is approached; modes must become progressively 
more anharmonic and smaller excitation energies should 
suffice to force the system into new ground states. Not 
only does the jamming/unjamming transition create crit- 
ical behavior in the harmonic properties of the solid [6] 
but also in its anharmonic response. The stress and tem- 
perature domain over which harmonic theory applies and 
anharmonic effects, such as those studied here, can be ne- 
glected is controlled by the proximity to this transition. 

* * * 

This work was supported by DE-FG02-05ER46199 (A. 
L., N. X., and V. V.), DE-FG02-03ER46088 (S. N. and N. 
X.), NSF-DMR05-47230 (V. V.), and NSF-MRSEC DMR- 
0820054 (S. N.). 



REFERENCES 

[1] KiTTEL, C, Solid State Physics (Wiley Publishing Inc., 

New York) 2004. 
[2] SiLBERT L. E., Liu A. J. and Nagel S. R., Phys. Rev. 

Lett, 95 (2005) 098301. 
[3] Zeravcic Z., Xu N., Liu A. J., Nagel S. R. and van 

Saarloos W., Europhys. Lett, 87 (2009) 26001. 
[4] SoMFAi, E., VAN Hecke, M., Ellendroek, W. G., 

Shundyak, K. and van Saarloos, W., Phys. Rev. E, 

75 (2007) 020301. 
[5] SiLBERT L. E., Liu A. J. and Nagel S. R., Phys. Rev. E, 

79 (021308) 2009. 



[6] O'Hern C. S., SiLBERT L. E., Liu A. J. and Nagel S. 

R., Phys. Rev. E, 68 (2003) 011306. 
[7] Xu N., ViTELLi v., Wyart M., Liu A. J. and Nagel S. 

R., Phys. Rev. Lett, 102 (2009) 038001. 
[8] ViTELLi v., Xu N., Wyart M., Liu A. J. and Nagel S. 

R., arXiv:0908.2176 (2009). 
[9] Press W. H., Flannery B. P. Teukolsky S. A. and 

Vetterling W. T., Numerical Recipes in FORTRAN, 

Vol. 77 (Cambridge University, New York) 1986. 
[10] Wyart M., Nagel S. R. and Witten T. A., Europhys. 

Lett, 72 (2005) 486; Wyart M., Silbert L. E., Nagel 

S. R. and Witten T. A., Phys. Rev. E, 72 (2005) 051306; 

Wyart M., Annales de Phys., 30(3) (2005) 1. 
[11] Xu N., Wyart M., Liu A. J. and Nagel S. R., Phys. 

Rev. Lett, 98 (2007) 175502. 
[12] http:/ /www. caam.rice.edu/software/ARPACK. 
[13] Xu N., Frenkel,D., et al, to be published (2009). 
[14] BuCHENAU U., et al, Phys. Rev. B, 46 (1992) 2798. 
[15] SCHOBER H. R. and Oligschleger C, Phys. Rev. B, 53 

(1996) 11469. 

[16] Schober H. R. and Ruocco G., Phil. Mag., 84 (2004) 
1361. 

[17] Fabian J. and Allen P. B., Phys. Rev. Lett, 79 (1997) 

1885; Phys. Rev. Lett, 77 (1996) 3839. 
[18] Biswas R., et al, Phys. Rev. Lett, 60 (1988) 2280. 
[19] Taraskin S. N. and Elliott S. R., Europhys. Lett, 39 

(1997) 37; Phys. Rev. B, 59 (1999) 8572. 

[20] WIDMER-COOPER A., et al. Nature Phys., 4 (2008) 711. 
[21] Zeravcic Z., van Saarloos W. and Nelson D. R., 

Europhys. Lett, 83 (2008) 44001. 
[22] BuCHENAU U., et al, Phys. Rev. Lett, 77 (1996) 659; 

Vainer Yu. G., et al., Phys. Rev. Lett, 97 (2006) 

185501. 

[23] Falk M. L. and Langer J. S., Phys. Rev. E, 57 (1998) 
7192; Langer J. S., Phys. Rev. E, 77 (2008) 021502. 



